library("effects")
library("foreign")
library("gmodels")
library("gplots")
library("nlme")
library("nnet")
library("stargazer")
library(MASS)
library(rockchalk)
library(pscl)
library(car)
library(DescTools)
library(psych)

d1 = read.dta("genpop replication.dta")
d1$race=relevel(d1$race, ref="White")
d1$dideo=relevel(d1$dideo, ref="Moderate")
d1$PartyID=relevel(d1$PartyID, ref="Independent")
d1$educ=relevel(d1$educ, ref="Some College")
d1$Age=Recode(d1$Age, "-936='18'")
d1$Age=as.numeric(d1$Age)
d1$Age=d1$Age+17
##Logit Models
mod.40 = multinom(Q23_1~Extra+Open+consc+Agree+emotstab+educ+income+race+dideo+PartyID+gender+Age, data=d1)
mod.50 = polr(Q23_1~Extra+Open+consc+Agree+emotstab+educ+income+race+dideo+PartyID+gender+Age, data=d1)
stargazer(mod.40,mod.50, type='html', style='ajps')
##Reviewer 2 would like to see the models without these controls
mod.140 = multinom(Q23_1~Extra+Open+consc+Agree+emotstab, data=d1)
mod.150 = polr(Q23_1~Extra+Open+consc+Agree+emotstab, data=d1)
stargazer(mod.140,mod.150, type='html', style='ajps')
##Calculate Pseudo R-squared
##Models
PseudoR2(mod.150, "all")
mod.400 = mnp(Q23_1~Extra+Open+consc+Agree+emotstab+educ+income+race+dideo+PartyID+gender+Age, data=d1)
mod.500 = polr(Q23_1~Extra+Open+consc+Agree+emotstab+educ+income+race+dideo+PartyID+gender+Age, data=d1)
mod.41 = multinom(Q23_1~educ+income+race+gender+Age+Extra+Open+Agree+consc+emotstab+infpp+polcar+businessC+socialC+CivicDuty+Candidate+lovePols+issue+supprty+servePub+Age, data=d1)
mod.51=polr(Q23_1~educ+income+race+gender+Age+Extra+Open+Agree+consc+emotstab+infpp+polcar+businessC+socialC+CivicDuty+Candidate+lovePols+issue+supprty+servePub+Age, data=d1)
stargazer(mod.40,mod.41, type='html', style='ajps')
##Grab Effects for Plotting Figure 1
effagree = as.data.frame(effect(c('Agree'), mod=mod.51))
effcons = as.data.frame(effect(c('consc'), mod=mod.51))
effextra = as.data.frame(effect(c('Extra'), mod=mod.51))
effopen = as.data.frame(effect(c('Open'), mod=mod.51))
effneuro = as.data.frame(effect(c('emotstab'), mod=mod.51))
write.csv(effagree, "Predicted Probability (agree).csv")
write.csv(effcons, "Predicted Probability (cons).csv")
write.csv(effextra, "Predicted Probability (extra).csv")
write.csv(effopen, "Predicted Probability (open).csv")
write.csv(effneuro, "Predicted Probability (emotstab).csv")
##Create Kernel Density Plots for Each of the Big Five Traits in both surveys and for those interested in running


d2 = read.dta("AMOS replication.dta")
d1$q23=as.numeric(d1$Q23_1)
d1rh=subset(d1,q23>1 )
d2rh=subset(d2, progamb_runhigher>1)
d1male=subset(d1, gender==c("Male"))
d1female=subset(d1, gender==c("Female"))
d2$gender = factor(d2$femaleprop, labels=c("Male", "Female"))
d2male=subset(d2, gender==c("Male"))
d2female=subset(d2, gender==c("Female"))
##Make the Agreeableness plot
plot(d2ad, col=3,lty=3, main="Agreeableness",ylim=c(0,2), xlim=c(0,3.5),lwd=2)
lines(d1rhad, col=2,lty=2, lwd=2)
lines(d1ad, lty=1)
#lines(d2rhad, col=4,lty=4, lwd=2)
legend("topleft", legend=c("Black = General Population", "Red = General Population Interested in Elected Office", "Green = Municipal Officials"))
##Grab the Densities for Extraversion
d1ed=density(d1$Extra, adjust=2, na.rm=T)
d1rhed=density(d1rh$Extra,adjust=2, na.rm=T)
d2ed=density(d2$extra4pt,adjust=2, na.rm=T)
d2rhed=density(d2rh$extra4pt,adjust=2, na.rm=T)
d1med=density(d1male$Extra,adjust=2, na.rm=T)
d1fed=density(d1female$Extra,adjust=2, na.rm=T)
d2med=density(d2male$extra4pt,adjust=2, na.rm=T)
d2fed=density(d2female$extra4pt,adjust=2, na.rm=T)

##Make the Extraversion plot
plot(d1ed, main="Extraversion", lty=1, ylim=c(0,2), xlim=c(0,3.5))
lines(d1rhed, col=2,lty=2, lwd=2)
lines(d2ed, col=3,lty=3,main="Extraversion", lwd=2)
#lines(d2rhed, col=4,lty=4, lwd=2)
legend("topleft", legend=c("Black = General Population", "Red = General Population Interested in Elected Office", "Green = Municipal Officials"))
##Grab the Densities for Conscientiousness
d1cd=density(d1$consc,adjust=2, na.rm=T)
d1rhcd=density(d1rh$consc,adjust=2, na.rm=T)
d2cd=density(d2$consc4pt,adjust=2, na.rm=T)
d2rhcd=density(d2rh$consc4pt,adjust=2, na.rm=T)
d1mcd=density(d1male$consc,adjust=2, na.rm=T)
d1fcd=density(d1female$consc,adjust=2, na.rm=T)
d2mcd=density(d2male$consc4pt,adjust=2, na.rm=T)
d2fcd=density(d2female$consc4pt,adjust=2, na.rm=T)

##Make the Conscientiousness plot
plot(d2cd, col=3,lty=3,main="Conscientiousness",ylim=c(0,2), xlim=c(0,3.5), lwd=2)
lines(d1rhcd, col=2,lty=2, lwd=2)
lines(d1cd, main="Conscientiousness", lty=1)
#lines(d2rhcd, col=4,lty=4, lwd=2)
legend("topleft", legend=c("Black = General Population", "Red = General Population Interested in Elected Office", "Green = Municipal Officials"))
##Grab the Densities for Openness
d1od=density(d1$Open,adjust=2, na.rm=T)
d1rhod=density(d1rh$Open,adjust=2, na.rm=T)
d2od=density(d2$open4pt,adjust=2, na.rm=T)
d2rhod=density(d2rh$open4pt,adjust=2, na.rm=T)
d1mod=density(d1male$Open,adjust=2, na.rm=T)
d1fod=density(d1female$Open,adjust=2, na.rm=T)
d2mod=density(d2male$open4pt,adjust=2, na.rm=T)
d2fod=density(d2female$open4pt,adjust=2, na.rm=T)

##Make the Openness plot
plot(d2od, col=3,lty=3,main="Openness",ylim=c(0,2), xlim=c(0,3.5), lwd=2)
lines(d1rhod, col=2,lty=2, lwd=2)
lines(d1od,main="Openness", lty=1)
#lines(d2rhod, col=4,lty=4, lwd=2)
legend("topleft", legend=c("Black = General Population", "Red = General Population Interested in Elected Office", "Green = Municipal Officials"))
##Grab the Densities for Emotional Stability
d1od=density(d1$emotstab,adjust=2, na.rm=T)
d1rhod=density(d1rh$emotstab,adjust=2, na.rm=T)
d2od=density(d2$stable4pt,adjust=2, na.rm=T)
d2rhod=density(d2rh$stable4pt,adjust=2, na.rm=T)
d1mod=density(d1male$emotstab,adjust=2, na.rm=T)
d1fod=density(d1female$emotstab,adjust=2, na.rm=T)
d2mod=density(d2male$stable4pt,adjust=2, na.rm=T)
d2fod=density(d2female$stable4pt,adjust=2, na.rm=T)
##Make the Emotional Stability plot
plot(d2od, col=3,lty=3,main="Emotional Stability",ylim=c(0,2), xlim=c(0,3.5), lwd=2)
lines(d1rhod, col=2,lty=2, lwd=2)
lines(d1od,main="Emotional Stability", lty=1)
#lines(d2rhod, col=4,lty=4, lwd=2)
legend("topleft", legend=c("Black = General Population", "Red = General Population Interested in Elected Office", "Green = Municipal Officials"))

##AMOS Models
d2$runhigher = factor(d2$progamb_runhigher)
d2$progamb=as.numeric(d2$runhigher)
d2$progamb_winlegis_1[d2$progamb_winlegis_1==-99] <- NA
mod.51 = multinom(progamb_runhigher~extra+ consc+ stable+ agree+ open, data=d2)
mod.52 = multinom(progamb_runhigher~extra+ consc+ stable+ agree+ open+termlimits+ closevote+ partisanelect+ progamb_current+ tenure+ progamb_similar_1+ progamb_winlegis_1, data=d2)
mod.54 = multinom(progamb_runhigher~extra+ consc+ stable+ agree+ open+termlimits+ closevote+ partisanelect +progamb_current+ tenure+ progamb_similar_1+ progamb_winlegis_1*agree, data=d2)
stargazer(mod.51,mod.52,mod.54, type='html')
##Ordered logit
mod.50 = polr(runhigher~extra+ consc+ stable+ agree+ open+gender+termlimits+ closevote+ partisanelect +progamb_current+ tenure+ progamb_similar_1+ progamb_winlegis_1, data=d2)
PseudoR2(mod.50, "all")
mod.55 = polr(runhigher~extra+ consc+ stable+ agree+ open+gender+termlimits+ closevote+ partisanelect +progamb_current+ tenure+ progamb_similar_1+ progamb_winlegis_1*agree, data=d2)
mod.55.ef = as.data.frame(effect(c('agree:progamb_winlegis_1'),mod=mod.55))
write.csv(mod.55.ef, "mod.55.ef.csv")
write.csv(mod.55.ef, "Figure3(1).csv")
PseudoR2(mod.55, "all")
##Reviewer 4 would like to see the models on a 4 point scale
mod.71=polr(runhigher~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt, data=d2)
mod.72 = polr(runhigher~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt+ closevote+ tenure+ progamb_current+termlimits+ partisanelect+ tenure+ progamb_similar_1+ progamb_winlegis_1, data=d2)
mod.73 = polr(runhigher~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt+ closevote+ tenure+ progamb_current+termlimits+ partisanelect+ tenure+ progamb_similar_1+ progamb_winlegis_1+gender, data=d2)
stargazer(mod.71,mod.72,mod.73, type='html')
mod.72.o.ef = as.data.frame(effect(c('open4pt'),mod=mod.72))
mod.72.c.ef = as.data.frame(effect(c('consc4pt'),mod=mod.72))
mod.72.e.ef = as.data.frame(effect(c('extra4pt'),mod=mod.72))
mod.72.a.ef = as.data.frame(effect(c('agree4pt'),mod=mod.72))
mod.72.n.ef = as.data.frame(effect(c('stable4pt'),mod=mod.72))
write.csv(mod.72.o.ef, "mod.72.o.ef.csv")
write.csv(mod.72.c.ef, "mod.72.c.ef.csv")
write.csv(mod.72.e.ef, "mod.72.e.ef.csv")
write.csv(mod.72.a.ef, "mod.72.a.ef.csv")
write.csv(mod.72.n.ef, "mod.72.n.ef.csv")

##Attractiveness of Higher Office
mod.100 = polr(office~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt, data=d2)
mod.101 = polr(office~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt+ closevote+ tenure+ progamb_current+termlimits+ partisanelect+ tenure+ progamb_similar_1+ progamb_winlegis_1+gender, data=d2)
mod.111 = multinom(office~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt+ closevote+ tenure+ progamb_current+termlimits+ partisanelect+ tenure+ progamb_similar_1+ progamb_winlegis_1+gender, data=d2)
stargazer(mod.100,mod.101, type='html', style='ajps')

##Interaction
##Reviewer 1 If the probability of winning question asks specifically about the probability of winning a state legislative seat, why is the interaction with agreeableness done with the generic DV from Table 2 and not the specific question about running for state office as in Table 3
mod.112 = runhigher(office~extra+consc+stable+agree+open+closevote+tenure+progamb_current+termlimits+partisanelect+progamb_similar_1+progamb_winlegis_1*agree, data=d2)
stargazer(mod.112, type='text', style='apsr')
plot(allEffects(mod.112), select="agree:progamb_winlegis_1", xlab="Probability of Winning", ylab="Predicted Probability", main="")
##Ordered Logit
mod.113 = polr(runhigher~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt+ closevote+ tenure+ progamb_current+termlimits+ partisanelect+ tenure+ progamb_similar_1+ progamb_winlegis_1+gender, data=d2)
mod.114 = polr(runhigher~extra4pt+ open4pt+agree4pt+ consc4pt+ stable4pt+ closevote+ tenure+ progamb_current+termlimits+ partisanelect+ tenure+ progamb_similar_1+ +gender+progamb_winlegis_1*agree4pt, data=d2)
PseudoR2(mod.113, "all")
PseudoR2(mod.114, "all")
allEffects(mod.114)
allEffects(mod.113)
effect('agree4pt:progamb_winlegis_1', mod.114, xlevels=list(agree4pt=c(0,2,2.75,3)))
plot(effect('agree4pt:progamb_winlegis_1', mod.114, xlevels=list(agree4pt=c(0,2,2.75,3))))

